**************************************************
* This code file set up the data and produces the
* figures and tables for the article:
*
* "THE QUALITY AND EFFICIENCY OF PUBLIC AND 
* PRIVATE FIRMS: EVIDENCE FROM AMBULANCE SERVICES"
*
* Accepted for publication at the
* Quarterly Journal of Economics
*
* By Daniel Knutsson and Bj�rn Tyrefors
* Research Institute of Industrial Economics
* Started 2018, final version 2022-02-10.
**************************************************


**************************************************
* Note. Complementary data is available in code file
* statistics_sweden.do provided in the replication package.
* This file shows how we have used matched
* employer-employee data that we use
* in the paper.
* ************************************************

**************************************************
* All data used in this code file has been provided by 
* SOS Alarm AB (coordinates) and the
* VAL database maintained by the County of Stockholm.
**************************************************
* Data sets (source):
*
* ambulance.dta (VAL)
* hospital_care.dta (VAL)
* primary_care.dta (VAL)
* flytt.dta (VAL)
* ambu_avl_alla.xlsx (VAL)
* SLL_InsatsNrMedPositionStatus.txt (SOSAlarm AB)
* SLL_InsatsNrMedPosition.txt (SOSAlarm AB)
**************************************************

* Programs required:
ssc install ftools, replace
ssc install geodist, replace
ssc install geo2xy, replace
ssc install spmap, replace
ssc install tmap, replace
ssc install shp2dta, replace
ssc install mif2dta, replace
ssc install reghdfe, replace
ssc install esttab, replace
ssc install ritest, replace
ssc install psacalc, replace

****************
* BEGIN CODE
****************

* Several datasets had to be recoded to unicode

cd C:\Data

unicode analyze ambulans.dta
unicode encoding set Windows-1252
unicode translate ambulans.dta

clear

cd C:\Data

unicode analyze hospital_care.dta
unicode encoding set Windows-1252
unicode translate hospital_care.dta

* Compress large file
use hospital_care.dta, clear
compress
save hospital_care.dta, replace

clear

*****************************************
* Start processing of data
*****************************************
* Hospital care visits 2007-2008
* Pre-treatment health variables

* 1. Number of hospital visits during these two years
* 2. Collect most frequent icd diagnosis for diseases in data
* 3. From these, create variables indicating a diagnosis 2007-2008

* Make the dataset smaller (hospital visits data)

cd C:\Data
use hospital_care.dta, clear

* Year variable
gen ar=substr(indat,1,4)
keep if ar=="2007" | ar=="2008"

* Keep identifier and diagnoses codes (at most 12 diagnoses)
keep NyID DIAG*

* Save this much smaller file for future use
save hc0708.dta, replace


* Start hospital visit data

use hc0708.dta, clear

* Number of hospital visits 2007-2008

gen one=1
egen ph_vis=sum(one), by(NyID)

* Take out most frequent diagnosis
* Make dataset long on the 12 diagnose variables

* Identify unique observation
seq id, f(1)

* Reshape 
reshape long DIAG, i(id) j(Nr)

* Drop if no diagnosis (i.e. if fewer than 12 diagnoses recorded)
drop if DIAG==""

* Drop a few duplicates
egen tag=tag(DIAG NyID)
keep if tag==1

* Most common diagnoses in data
egen sum=sum(one), by(DIAG)

* Here, I have manually checked the 22 most common

gen dia=.
gen hyp=.
gen heart=.
gen kol=.
gen angi=.
gen opi=.
gen dep=.
gen alc=.
gen reu=.
gen pain=.

gen anemi=.
gen flim=.
gen heart2=.
gen uvi=.
gen hyty=.
gen hylip=.
gen stroke=.
gen pneu=.

gen pros=.
gen back=.
gen anx=.
gen pros2=.

gen new=""

replace new=substr(DIAG,1,3)

* Define variables based on the most common diagnoses (ICD 10 codes, first three digits)

replace dia=1 if new=="E10" | new=="E11" | new=="E12" | new=="E13" | new=="E14" 
replace hyp=1 if new=="I10" | new=="I11" | new=="I12" | new=="I13" | new=="I14" | new=="I15" 
replace heart=1 if new=="I50"
replace kol=1 if new=="J40" | new=="J41" | new=="J42" | new=="J43" | new=="J44" | new=="J45" | new=="J46" | new=="J47"
replace angi=1 if new=="I20"
replace opi=1 if new=="F11"
replace dep=1 if new=="F32"
replace alc=1 if new=="F10"
replace reu=1 if new=="M05"
replace pain=1 if new=="R52"

replace anemi=1 if new=="D64"
replace flim=1 if new=="I48"
replace heart2=1 if new=="I25"
replace uvi=1 if new=="N39"
replace hyty=1 if new=="E03"
replace hylip=1 if new=="E78"
replace stroke=1 if new=="I63"
replace pneu=1 if new=="J18"

replace pros=1 if new=="C61"
replace pros2=1 if new=="N40"
replace back=1 if new=="M54"
replace anx=1 if new=="F41"

lab var heart "Heart Failure"
lab var heart2 "Chronic Heart Disease"

foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace `x'=0 if `x'==.
}

* Take the data back to the individual level
collapse (sum) dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 (mean) ph_vis , by(NyID)

* Make the variable an indicator variable instead of number of times a diagnosis was recorded
foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace `x'=1 if `x'>0
}


lab var ph_vis "Hospital Visits 2007-2008"

save phvis.dta, replace


************************************************
* Non-hopspital care (primary care visits) from the VAL databse
* Do the same for this data
************************************************

* Load data in pieces and translate to unicode format
* Data contains long string variables and many observations

cd C:\Data

set more off

local y=1

local k=1

forval x=10000000(10000000)99452826 {

use lkomm basomr bdat btyp DIAG1 DIAG2 DIAG3 DIAG4 DIAG5 ATG1 ATG2 ATG3 ATG4 ATG5 typ lkf NyID mosaicgrupp ar in `k'/`x' using primary_care.dta, clear

compress

save pc`y'.dta, replace

clear

unicode encoding set Windows-1252

unicode translate pc`y'.dta

local k=`x'+1

local y=`y'+1

}

***************************************************
* Primary care visits in 2007-2008
* pc8.dta and pc9.dta cover 2007-2008
***************************************************

cd C:\Data

use pc8.dta, clear

append using pc9.dta

destring ar, replace

drop if ar==2006 | ar==2009

* Code is similar to above using hospital visits, see comments there

keep NyID DIAG* 

seq id, f(1)

gen one=1

egen pc_vis=sum(one), by(NyID)

reshape long DIAG, i(id) j(Nr)

drop if DIAG==""

egen tag=tag(DIAG NyID)

keep if tag==1

* egen sum=sum(one), by(DIAG)

gen dia=.
gen hyp=.
gen heart=.
gen kol=.
gen angi=.
gen opi=.
gen dep=.
gen alc=.
gen reu=.
gen pain=.

gen anemi=.
gen flim=.
gen heart2=.
gen uvi=.
gen hyty=.
gen hylip=.
gen stroke=.
gen pneu=.

gen pros=.
gen back=.
gen anx=.
gen pros2=.

gen new=""

replace new=substr(DIAG,1,3)
replace dia=1 if new=="E10" | new=="E11" | new=="E12" | new=="E13" | new=="E14" 
replace hyp=1 if new=="I10" | new=="I11" | new=="I12" | new=="I13" | new=="I14" | new=="I15" 
replace heart=1 if new=="I50"
replace kol=1 if new=="J40" | new=="J41" | new=="J42" | new=="J43" | new=="J44" | new=="J45" | new=="J46" | new=="J47"
replace angi=1 if new=="I20"
replace opi=1 if new=="F11"
replace dep=1 if new=="F32"
replace alc=1 if new=="F10"
replace reu=1 if new=="M05"
replace pain=1 if new=="R52"

replace anemi=1 if new=="D64"
replace flim=1 if new=="I48"
replace heart2=1 if new=="I25"
replace uvi=1 if new=="N39"
replace hyty=1 if new=="E03"
replace hylip=1 if new=="E78"
replace stroke=1 if new=="I63"
replace pneu=1 if new=="J18"

replace pros=1 if new=="C61"
replace pros2=1 if new=="N40"
replace back=1 if new=="M54"
replace anx=1 if new=="F41"


foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace `x'=0 if `x'==.
}

collapse (sum) dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 (mean) pc_vis, by(NyID)

foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace `x'=1 if `x'>0
rename `x' VC`x'
}


lab var pc_vis "Primary Care Visits 2007-2008"

save pc_0708.dta, replace


**********************************************
* Mortality data (the VAL database)
**********************************************
cd C:\Data

use flytt.dta, clear

* Code for dead
keep if hkod=="D"

drop hkod

* Create death date variable
rename hdat ddate

* Drop a few duplicates
duplicates t NyID, gen(tag)
tab tag
drop if tag>0

compress

save death.dta, replace


*********************************************
* Mortality data no 2, complemented in 2021
* This data was complemented from the VAL database
* to assure that we cover 3-year mortality
* up until 2016 (mortality data until 2019)
*********************************************

cd C:\Data

import excel using "ambu_avl_alla.xlsx", first clear

drop Nyhdat3

* No duplicates

save "nymort.dta", replace

************************************************************************
* Coordinates provided by SOSAlarm AB (ambulance dispatch firm)
************************************************************************

* First GPS coordinates 

cd C:\Data

import delimited using SLL_InsatsNrMedPosition.txt, varnames(1) delim("|")

drop datum

save ambco.dta, replace

clear


* Also save projected coordinates (x,y)

cd C:\Data

import delimited using SLL_InsatsNrMedPosition.txt, delim("|") varnames(1)

drop datum

* Coordinates:
* Clean GPS coordinates to conform with Stata

foreach z in x y {
gen first=substr(patient_`z',1,2)
gen second=substr(patient_`z',4,15)

gen p`z'k=first+"."+second
drop first second
}

* If missing string
replace pxk="" if pxk=="."
replace pyk="" if pyk=="."

* Check for string/numeric
gen byte notnumeric = real(pxk)==.

* Take away commas
replace pxk = subinstr(pxk, "," , "" ,.) 
replace pyk = subinstr(pyk, "," , "" ,.) 

* Replace to comma
gen lat=subinstr(pxk,".",",",.)
gen lon=subinstr(pyk,".",",",.)

destring pxk, gen(px) float
destring pyk, gen(py) float

format px %20.0g
format py %20.0g

replace px=. if px==0
replace py=. if py==0

* Project (x,y) GCS_WGS_1984
geo2xy py px , gen(pty ptx)  project(mercator_sphere)

format ptx %20.0g
format pty %20.0g

keep insatsnr pty ptx py px

* Restrict inclusion to values within a square
* covering Stockholm County
* 3333 + 118 deleted with out of range coordinates

drop if py<17.14 | py>19.158
drop if px<58.8 | px>60.31

save corsamp.dta, replace


**********************************************
* Time stamps from SOSAlarm AB
* Additional ambulance responeses (during the transport)
* that were not loaded to the VAL database but saved
* at SOSAlarm AB. We complemented these data to have better coverage.
* Can be merged by unique dispatch number. 
**********************************************

cd C:\Data

import delimited using SLL_InsatsNrMedPositionStatus.txt, varnames(1) delim("|")

save ambst.dta, replace

clear


**********************************************
* Hospital visits, admission and diagnoses
**********************************************
cd C:\Data

use "hospital_care.dta", clear
gen adate=date(indat,"YMD")
gen udate=date(utdat,"YMD")

sort NyID adate udate
drop if udate==adate[_n-1] & NyID==NyID[_n-1] & adate==adate[_n-1]
egen tag=tag(NyID adate)
keep if tag==1

gen inl=1

keep adate NyID inl DIAG*

save "hosp2.dta", replace


************************************************************************
* Start setting up data on Ambulances
* Main ambulance data file from VAL database
************************************************************************

cd C:\Data

use ambulans.dta, clear

* Date variables
gen date=date(manad, "YM")
gen month=month(date)
gen year=year(date)

* Merge with additional datasets

* Patient coordinates (projected)
merge 1:1 insatsnr using corsamp.dta
drop if _merge==2
drop _merge

* New status, time-stamps, from SOSAlarm AB
merge 1:1 insatsnr using ambst.dta
drop if _merge==2
drop _merge

* Mortality data

merge m:1 NyID using death.dta
drop if _merge==2
drop _merge

* Primary care visits 2007-2008

merge m:1 NyID using pc_0708.dta
replace pc_vis=0 if _merge==1
drop if _merge==2
drop _merge

* Hospital care visits 2007-2008

merge m:1 NyID using phvis.dta
replace ph_vis=0 if _merge==1
drop if _merge==2
drop _merge


************************************************
* Helicopter variables to drop
************************************************

drop flygtid flygframpat flygstart hkpframtid

****************************************************
* Generate variables
****************************************************

****************************************************
* Multiple ambulances. 3 ways of defining a multiple
* ambulance dispatch.
* Several patients or very severe condition
****************************************************

****************************************************
* 1. Multiple patients
****************************************************

* Ambulance dispatch number, several patients or
* potentially several patients.
* Equal numbers before ":" indicate same event but 
* different patient.

gen one=1
gen numb1=substr(insatsnr,1,4)
gen numb2=substr(insatsnr,5,.)
gen numb3=substr(numb2,1,strpos(numb2,":")-1)
gen numb=numb1+numb3

* Sum number of units per same event
egen nugr=group(numb)
egen nusu=sum(one), by(nugr)
tab nusu
drop numb1 numb2 numb3 nugr numb

****************************************************
* 2. Multiple units per patient
****************************************************
* Indicator for multiple ambulance dispatch
* If more ambulances to same patient, only the
* first ambulance gets a dispatch number.
* The other units are recorded on the same observation.
* I use these variables to define more than one unit 
* dispatched.

gen more_amb=0
foreach x in ANKTIDENH2 AVFARDENH2 FRAMENH2 KVITTENSENH2 MOBITEXENH2 ENHET2 LEV2 {
replace more_amb=1 if `x'!=""
}

****************************************************
* 3. Multiple units per patient
****************************************************
* Emergency car assisted. Precoded indicator if 
* specialist unit accompanied the ambulance.

gen akut=0
replace akut=1 if AKUTBILMED1=="J"


****************************************************
* Variable creation
****************************************************

* Date of transport
gen adate=date(datum, "YMD")

* Female
gen female=1 if koen=="K"
replace female=0 if koen=="M"


* Ambulance area for pickup (h�mtomr�de)
* Each unit has a local area designated to its unit code
* Here, I clean that data to comply with unit codes.
* Mainly by changing hopital descriptions to unit codes.

gen first=substr(hamtomr,1,4)
gen second=substr(first,2,3)
gen new=upper(substr(first,1,1))
gen apar=new+second

replace apar="DS" if apar=="DS." | apar=="DSB"  | apar=="Ds" | apar=="Ds."
replace apar="KS" if apar=="KS." | apar=="KSB"  | apar=="KSN" | apar=="KSR" | apar=="KST" | apar=="KSn" | apar=="KSt" | apar=="Ks" | apar=="Ks." | apar=="Ksb" | apar=="Ksn" | apar=="Ksr" | apar=="Kst"
replace apar="NTS" if apar=="Nts" 
replace apar="SGS" if apar=="SGS." | apar=="Sgs"  | apar=="Sgs."
replace apar="S�S" if apar=="S�s" | apar=="Z�S"
replace apar="HS" if apar=="HSB" | apar=="Hs" | apar=="Hsb"
replace apar="STS" if apar=="Sts"

replace apar="A921" if apar=="KS"
replace apar="A907" if apar=="SGS"
replace apar="B911" if apar=="NTS"
replace apar="B931" if apar=="STS"
replace apar="A901" if apar=="S�S"
replace apar="A981" if apar=="HS"
replace apar="A971" if apar=="DS"

replace apar="BMA" if apar=="BMa" | apar=="Bma" | apar=="ABRO"
replace apar="DAL" if apar=="ADAL" 
replace apar="SAB" if apar=="ASAB"  | apar=="SBS"  | apar=="Sbs"

replace apar="A909" if apar=="BMA"
replace apar="A903" if apar=="DAL"
replace apar="A907" if apar=="SAB"

* Group variable for this variable
egen gr=group(apar)

* If fewer than 5000 patients in a code, drop.
* These are descriptions of the patient that 
* could not be assigned an ambulance (district) code

egen sum=sum(one), by(gr)
replace apar="" if sum<5000

drop gr sum one

* Variable "apar" is the ambulance area of patient


* Treatment measure
gen private=1
replace private=0 if LEV1=="AISAB"

* Priority to patient and to hospital
destring priout, replace
destring prioin, replace

* Patients that were not transported to hospital (138 000 obs)

gen pmed=1
replace pmed=0 if patmed=="J"
replace pmed=. if patmed==""


* Ambulance last year (if any)
sort NyID adate
gen ambly=0
replace ambly=1 if NyID==NyID[_n-1] & adate<=adate[_n-1]+365


********************************************
* Health variables 2007-2008
********************************************
* For these variable, I include both in-patient and out-patient data
* on diagnoses. If any diagnosis, indicator is one.

foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace `x'=1 if VC`x'==1
replace `x'=0 if `x'==.
}

* Trunkate out-patient visits "pc_vis" to 400
* There are extreme outliers in this variable.
* Some patients may have had visits daily, or even more.

replace pc_vis=200 if pc_vis>200

* Any in-hopital or out-hospital visits variables
gen vcvi=pc_vis
replace vcvi=1 if vcvi>1

gen hvis=ph_vis
replace hvis=1 if hvis>1


* Precoded variable if the patient resides at
* a nursing home.
gen sb=0
replace sb=1 if vessb!=""

* The ambulance crew assessment of the severity
* of the patient condition.
gen cond=allvarlgrad
replace cond="" if allvarlgrad=="V"
destring cond, replace

********************************************
* Variables - Fixed effects - geographical
********************************************

* Grid FE based on patient coordinates
* Create grids for around k*k number of cells 
* All cells will not have coverage of observations
* All cells will not be included in Stockholm Region

* I have already restricted the coordinates to a square
* around Stockholm Region.

set more off 

forval z=200(50)1000 {

local k=`z'
qui sum ptx
local ymin=r(min)
local ymax=r(max)
local yrange=`ymax'-`ymin'
local yjump=`yrange'/`k'

qui sum pty
local xmin=r(min)
local xmax=r(max)
local xrange=`xmax'-`xmin'
local xjump=`xrange'/`k'

gen fe`k'=.
local s=1
forval x=`xmin'(`xjump')`xmax' {
local new=`x'+`xjump'
replace fe`k'=`s' if pty>=`x' & pty<`new'
local s=`s'+1
}

qui sum fe`k'
local grm=r(max)
local z=0
forval y=`ymin'(`yjump')`ymax' {
local new=`y'+`yjump'
local new2=`grm'*`z'
replace fe`k'=`new2'+fe`k' if ptx>=`y' & ptx<`new'
local z=`z'+1
}
}


*****************************************
* Hospital fixed effects
*****************************************

gen hosp=.
local k=1
replace hosp=`k' if avlstalle=="Astrid Lindgrens Barnsjukhus"
replace hosp=`k' if avlstalle=="Karolinska Solna"
replace hosp=`k' if avlstalle=="Karolinska neurokir"
replace hosp=`k' if avlstalle=="Karolinska soln"
replace hosp=`k' if avlstalle=="Karolinska solna"
replace hosp=`k' if avlstalle=="Karolinska thorax"
replace hosp=`k' if avlstalle=="Karolnska solna"
local k=`k'+1
replace hosp=`k' if avlstalle=="Danderyds sjukhus"
local k=`k'+1
replace hosp=`k' if avlstalle=="Huddinge sjukhus"
replace hosp=`k' if avlstalle=="Huddinge barnsjukhus"
local k=`k'+1
replace hosp=`k' if avlstalle=="Norrt�lje sjukhus"
local k=`k'+1
replace hosp=`k' if avlstalle=="St G�rans sjukhus"
local k=`k'+1
replace hosp=`k' if avlstalle=="S�dersjukhuset"
replace hosp=`k' if avlstalle=="Sachsska Barnsjukhus"
local k=`k'+1
replace hosp=`k' if avlstalle=="S�dert�lje sjukhus"


********************************************
* Outcomes - Time stamps
********************************************

gen double dispatch=clock(MOBITEXENH1,"DMYhms")
gen double confirm=clock(KVITTENSENH1,"DMYhms") 
gen double arrivep=clock(FRAMENH1,"DMYhms") 
gen double gohosp=clock(AVFARDENH1,"DMYhms") 
gen double arriveh=clock(ANKTIDENH1,"DMYhms") 


* Extra data complemented from SOSAlarm AB (same as some above but more complete)
foreach x in t s k {
replace status_`x'=substr(status_`x',1,strpos(status_`x',".")-1)
}

gen double dispatch2=clock(status_t,"YMDhms") 
gen double arriveh2=clock(status_s,"YMDhms")
gen double done=clock(status_k,"YMDhms")

replace dispatch2=dispatch if dispatch2==.
replace arriveh=arriveh2 if arriveh==.

* Variables are truncated for unreasonably long times.
* See Online Appendix Section F.2 describing why.

*************************************
* Response time
*************************************
gen resp=(confirm-dispatch2)/1000
replace resp=. if resp<0
* Not more than 30 minutes to confirm (cut at 30 min)
replace resp=. if resp>1800

*************************************
* Time To Patient 
*************************************
gen top=(arrivep-confirm)/1000
replace top=. if top<0

* Max time to patient (1h)
replace top=. if top>3600

*************************************
* Time At Patient
*************************************
gen ptime=(gohosp-arrivep)/1000
replace ptime=. if ptime<0

* Max time with patient (4h)
replace ptime=. if ptime>14400

*************************************
* Time To Hospital
*************************************
gen toh=(arriveh-gohosp)/1000
replace toh=. if toh<0

* Max time to hospital 2h
replace toh=. if toh>7200

*************************************
* Getting ambulance in order/break time/coffee
*************************************
gen fika=(done-arriveh)/1000
replace fika=. if fika<0

* Max time at hospital 2h
replace fika=. if fika>7200

*************************************
* Patient queue time
*************************************
gen double create=clock(SKAPENH1,"DMYhms")
gen pqt=(dispatch2-create)/1000
replace pqt=. if pqt<0


* Truncate waiting time to 10 h
replace pqt=36000 if pqt>36000

replace pqt=. if pqt==0

********************************************
* Outcomes - Mortality
********************************************

* Death from k years after ambulance service

gen deathd=date(ddate,"YMD")

gen d1=0
replace d1=1 if deathd<=adate+1095

gen d2=0
replace d2=1 if deathd<=adate+730

gen d3=0
replace d3=1 if deathd<=adate+365

gen d4=0
replace d4=1 if deathd<=adate+180

gen d5=0
replace d5=1 if deathd<=adate+90

gen d6=0
replace d6=1 if deathd<=adate+30

gen d7=0
replace d7=1 if deathd<=adate+7

gen d8=0
replace d8=1 if deathd<=adate+1



************************************************
* Sample restrictions
************************************************

* Out of county transport
drop if utomlan=="J"

* Drop unclear provider (1000 obs)
drop if LEV1=="Extra" 

* Observations that should not be in data
drop if year==2008 | year>2016

* No information, priority 0 does not exist
replace priout=. if priout==0

* Defining real ambulances (not special units)
gen amb=0
replace amb=1 if (ENHET1TYP=="Akutambulans" | ENHET1TYP=="Cert.Amb." | ENHET1TYP=="Reservambulans")

* More 1 unit dispatches (only 1 patient)
gen more=1
replace more=0 if nusu==1 & akut==0 & more_amb==0

* Akutambulansbilar placerade
gen aba=0
replace aba=1 if STATION1=="Huddinge" | STATION1=="Sollentuna"

gen extra=0
replace extra=1 if STATION1=="Extra" | ENHET1TYP=="Extraambulans"

* Time and date variables using dispatch services time of
* response

* Day variable

gen t1=substr(SKAPENH1,1,strpos(SKAPENH1,":")-1)

* Time variable
gen t2=substr(SKAPENH1,strpos(SKAPENH1,":")+1,.)
* Hour of the day
gen t3=substr(t2,1,strpos(t2,":")-1)

drop t2

destring t3, replace
rename t3 ctime

gen double calld=date(t1,"DMY")
gen dow=dow(calld)

gen week=week(calld)


* Patient picked up in a private area (area defined by a private ambulance code)
gen apriv=1
replace apriv=0 if apar=="A901" | apar=="A903" | apar=="A902" | apar=="A907" | apar=="A908" | apar=="A909" | apar=="A921" | apar=="A924" | apar=="A931"

* Public ambulance district codes:
* S�dermalm (A901, A903)
* City (A902, A907)
* Bromma (A908)
* H�sselby (A909) (V�llingby)
* Solna (A921)
* Eker� (A924)
* Liding� (A931)

* Hospital transfer, use pick up area descriptions (hospital specific acronyms)
drop first second new
gen first=substr(hamtomr,1,4)
gen second=substr(first,2,3)
gen new=upper(substr(first,1,1))
gen apar2=new+second

replace apar2="DS" if apar=="DS." | apar=="DSB"  | apar=="Ds" | apar=="Ds."
replace apar2="KS" if apar=="KS." | apar=="KSB"  | apar=="KSN" | apar=="KSR" | apar=="KST" | apar=="KSn" | apar=="KSt" | apar=="Ks" | apar=="Ks." | apar=="Ksb" | apar=="Ksn" | apar=="Ksr" | apar=="Kst"
replace apar2="NTS" if apar=="Nts" 
replace apar2="SGS" if apar=="SGS." | apar=="Sgs"  | apar=="Sgs."
replace apar2="S�S" if apar=="S�s" | apar=="Z�S"
replace apar2="HS" if apar=="HSB" | apar=="Hs" | apar=="Hsb"
replace apar2="STS" if apar=="Sts"

* Variable describing if patient was picked up from a hospital
gen fhosp=0
replace fhosp=1 if apar2=="KS" | apar2=="NTS" | apar2=="SGS" | apar2=="S�S" | apar2=="HS" | apar2=="STS" | apar2=="DS" 

*******************************************************
* Clean data to keep only variables to be used
*******************************************************

keep private deathd pc_vis ph_vis amb extra aba priout fe200 fe250 fe300 fe350 fe400 fe450 fe500 fe550 fe600 fe650 fe700 fe750 fe800 fe850 fe900 fe950 fe1000 year month priout week dow ctime NyID vcvi hvis pc_vis ph_vis dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 ambly alder female sb pqt resp top ptime toh fika cond prioin pmed d1 d1 d2 d3 d4 d5 d6 d7 d8 apriv hosp fhosp more avlstalle adate apar STATION1 BEDTILLST1 BEDTILLST2 BEDTILLST3 AKUTBILMED1 AKUTBILMED2 ENHET2 ENHET3 dispatch2

lab var fhosp "Transported from hospital"
lab var hosp "Destination Hospital"
lab var apriv "Private district"

* Normalize variables
qui sum alder
gen asd=(alder-r(mean))/r(sd)

lab var asd "Age normalized"

qui sum pc_vis
gen pvis=(pc_vis-r(mean))/r(sd)

lab var pvis "Primary care visits normalized"

qui sum ph_vis 
gen hovs=(ph_vis-r(mean))/r(sd)

lab var hovs "Hospital visits normalized"

* Priority 1 to hospital
gen p1in=0
replace p1in=1 if prioin==1
replace p1in=. if prioin==.

* Patient in queue time (set max 10 hours)
replace pqt=. if pqt>=36000

**********************************************
* Diagnoses in the ambulance
**********************************************

gen ca=0
forval x=1/1 {
replace ca=1 if BEDTILLST`x'=="Hj�rtstopp - asystoli - PEA" | BEDTILLST`x'=="Hj�rtstopp - med framg�ngsrik HLR"  | BEDTILLST`x'=="Hj�rtstopp-ventrikelflimmer"
}

gen ami=0
forval x=1/1 {
replace ami=1 if BEDTILLST`x'=="Hj�rtinfarkt"
}

gen ospec=0
forval x=1/1 {
replace ospec=1 if BEDTILLST`x'=="Allm�nt ospecificerat" | BEDTILLST`x'=="Allm�nt �vrigt"
}

gen bstro=0
forval x=1/1 {
replace bstro=1 if BEDTILLST`x'=="Hj�rnbl�dning/hj�rninfarkt"
}

**************************************************
* Hospital admissions by days after ambulance transport
**************************************************

gen adate_real=adate
gen inl_all=0
gen dia1=""

local j=0

forval x=0/30 {
replace adate=adate+`j'
qui merge m:1 NyID adate using "C:\Users\Min Dator\Dropbox\Ambulans\Analysis\hosp2.dta"
drop if _merge==2
replace inl_all=1 if inl==1
rename inl inl`x'
replace inl`x'=0 if inl`x'==.
replace dia1=DIAG1 if dia1==""
drop _merge
local j=1
}

drop adate

rename adate_real adate

lab var inl_all "Any hospital admission within 30 days"
lab var kol "COPD"
lab var angi "Angina pectoris"
lab var opi "Opioid disease"
lab var dep "Depressive disease"
lab var alc "Alcohol disease"
lab var reu "Reumatic"
lab var pain "Cronic pain"
lab var anemi "Anemia"
lab var flim "Cardiac fibrillation"
lab var heart2 "Chronic heart disease"
lab var uvi "Urinary infection"
lab var hyty "Hyperthyroidism"
lab var hylip "Hyperlipidaemia"
lab var stroke "Stroke"
lab var pneu "Pneumonia"
lab var pros "Benign prostate cancer"
lab var back "Back pain"
lab var anx "Anxiety"
lab var pros2 "Malign prostate cancer"
lab var heart "Heart failure"
lab var heart2 "Chronic ischemic heart disease"


******************************************
* Interactions for heterogeneity analyses
******************************************

* Private and priority 1 to patient
gen p1=0
replace p1=1 if priout==1
gen popri=private*p1

* Private and female
gen spri=female*private

* Above age 70 and private
gen old=0
replace old=1 if alder>70
gen alpri=private*old

* Heart disease diagnosis and private
gen heartpri=heart2*private

* Angina diagnosis and private
gen anpri=angi*private

* COPD diagnosis and private
gen resppri=private*kol

* Urinary tract infection and private
gen uvipri=pain*private

* Depression and private
gen addpri=private*dep

drop DIAG* dia1

********************************************
* Admissions, same day and one day after
* The ambulance transport and admission
* decision can often be on different days.
* Take out ICD codes for intake diagnosis
********************************************

forval x=1/9 {
gen dia`x'=""
}

* Same day
qui merge m:1 NyID adate using "C:\Data\hosp2.dta"
drop if _merge==2
drop _merge
forval x=1/9 {
rename DIAG`x' diag`x'
}

* Day after
replace adate=adate+1
qui merge m:1 NyID adate using "C:\Data\hosp2.dta"
drop if _merge==2
drop _merge
forval x=1/9 {
replace diag`x'=DIAG`x' if diag`x'==""
}

drop DIAG*
replace adate=adate-1
drop inl
gen inl=0
replace inl=1 if diag1!=""

lab var inl "Admission up to one day after ambulance service"

* Intake diagnoses:
* First icd code only (use first letter codes with more 
* than 8000 main diagnoses. 0.1% av data lost)

gen temp=substr(diag1,1,1)
foreach z in A C E F G I J K M N R S T Z {
gen `z'icd2=0
replace `z'icd2=1 if temp=="`z'"
}

* More units to same patient and high priority
* Indication of severe event

gen more2=0
replace more2=1 if ENHET2!="" & priout==1

****************************************************
* Make all tables communicate with diseases
* Ambulance diagnoses grouped
****************************************************

* Respiratory
gen bresp=0
replace bresp=1 if BEDTILLST1=="Andning ospecificerat" | BEDTILLST1=="Andning �vrigt" | BEDTILLST1=="Andningsbesv�r" | BEDTILLST1=="Andningsbesv�r m framm. kropp" | BEDTILLST1=="Andningsbesv�r m framme. kropp" | BEDTILLST1=="Andningsbesv�r m fr�mmande kropp" | BEDTILLST1=="Andningsbesv�r m pip astma" | BEDTILLST1=="Andningsstillest�nd" | BEDTILLST1=="Pseudokrupp" | BEDTILLST1=="Pneumothorax" | BEDTILLST1=="Lunginflammation" | BEDTILLST1=="Lungemfysem" | BEDTILLST1=="Kronisk obstruktiv lungsjukdom" | BEDTILLST1=="Kronisk bronkit" | BEDTILLST1=="Influensa" | BEDTILLST1=="Epiglottit"

replace bresp=. if BEDTILLST1==""

* CVD
gen bcvd=0
replace bcvd=1 if BEDTILLST1=="Aortaaneurysm" | BEDTILLST1=="Centrala br�stsm�rtor - k�rlkramp" | BEDTILLST1=="Cirkulation" | BEDTILLST1=="Cirkulation," | BEDTILLST1=="Cirkulation, ospecificerat" | BEDTILLST1=="Cirkulation, �vrigt" | BEDTILLST1=="Embolier tromboser - art�r" | BEDTILLST1=="Embolier tromboser - ven" | BEDTILLST1=="Hj�rnbl�dning/hj�rninfarkt" | BEDTILLST1=="Hj�rnbl�dning/hj�rninfarkt, svalg, ma.." | BEDTILLST1=="Hj�rtinfarkt" | BEDTILLST1=="Hj�rtstopp - asystoli - PEA" | BEDTILLST1=="Hj�rtstopp - med framg�ngsrik HLR" | BEDTILLST1=="Hj�rtstopp-ventrikelflimmer" | BEDTILLST1=="Hj�rtsvikt" | BEDTILLST1=="Hypertoni" | BEDTILLST1=="Hypotoni" | BEDTILLST1=="Lung�dem" | BEDTILLST1=="Subarachnoidalbl�dning" | BEDTILLST1=="Transitorisk ischemisk attack" | BEDTILLST1=="entrala br�stsm�rtor - k�rlkrampC" | BEDTILLST1=="Rytm- eller retledningshinder"
replace bcvd=. if BEDTILLST1==""

* GI
gen bsur=0
replace bsur=1 if BEDTILLST1=="mage-tarm" | BEDTILLST1=="Sjukdomar i bukspottsk�rteln" | BEDTILLST1=="Mags�r" | BEDTILLST1=="Kirurgi �vrigt" | BEDTILLST1=="Kirurgi ospecificerat" | BEDTILLST1=="Ileus" | BEDTILLST1=="Br�stkorgssm�rtor" | BEDTILLST1=="Br�ck" | BEDTILLST1=="Bl�dning fr�n mun, svalg, mage-tarm" | BEDTILLST1=="Bl�dning fr�n mun, mage-tarm"
replace bsur=. if BEDTILLST1==""
 
* Psychyatry/drugs/intoxication
gen bpsy=0
replace bpsy=1 if BEDTILLST1=="Drog-/l�kemedelsmissbruk" | BEDTILLST1=="F�rgiftning med alkohol" | BEDTILLST1=="F�rgiftning med alkohol och l�kemedel" | BEDTILLST1=="F�rgiftning med f�do�mne inkl svamp" | BEDTILLST1=="F�rgiftning med l�kemedel" | BEDTILLST1=="F�rgiftning med narkotika" | BEDTILLST1=="F�rgiftning med tobak" | BEDTILLST1=="F�rgiftning ospecificerat" | BEDTILLST1=="F�rgiftning petroleumprod inkl etylen.." | BEDTILLST1=="F�rgiftning �vrigt" | BEDTILLST1=="Hallucinationer" | BEDTILLST1=="Krisreaktion" | BEDTILLST1=="Psykiatri ospecificerat" | BEDTILLST1=="Psykiatri �vrigt" | BEDTILLST1=="Psykiska symtom och sjukdomar"
replace bpsy=. if BEDTILLST1==""

*****************************************************
* Pre-determined health (2007-2008)
* Classify in similar terms as above
*****************************************************

* CVD
gen pcvd=0
replace pcvd=1 if heart==1 | heart2==1 | angi==1 | flim==1 | stroke==1

* Respiratory
gen presp=0
replace presp=1 if pneu==1 | kol==1

* Psychiatric/intoxication
gen ppsy=0
replace ppsy=1 if opi==1 | dep==1 | alc==1 | anx==1

* Pain/cancer
gen psur=0
replace psur=1 if back==1 | pain==1 | anemi==1 | pros==1 | pros2==1

foreach x in pcvd presp ppsy psur {
gen pri`x'=`x'*private
}

* Do male*private interaction
gen female_real=female
replace female=0
replace female=1 if female_real==0
drop spri
gen spri=female*private

* Change to "no hospital visit" for interaction
gen hvis2=0 if hvis==1
replace hvis2=1 if hvis==0
gen prihvis=hvis2*private

* Multiple diagnoses 2007-2008 (more than 1)
* For interaction/heterogeneity
gen multis=0
foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 {
replace multis=multis+1 if `x'==1
}

replace multis=0 if multis==1
replace multis=1 if multis>1

gen primult=private*multis

lab var female "Male"
lab var multis "No hospital visits 2007-2008"
compress

lab var pmed "Patient not conveyed"

**********************************************
* Add new mortality data
* Data for 2018-2019. Updated 2021
**********************************************

merge m:1 NyID using "C:\Data\nymort.dta"
drop if _merge==2
drop _merge

gen deathd2=date(HDAT,"YMD")

gen xd0=0
replace xd0=1 if deathd2<=adate+1460
replace xd0=. if year==2016

gen xd1=0
replace xd1=1 if deathd2<=adate+1095

gen xd2=0
replace xd2=1 if deathd2<=adate+730

gen xd3=0
replace xd3=1 if deathd2<=adate+365

gen xd4=0
replace xd4=1 if deathd2<=adate+180

gen xd5=0
replace xd5=1 if deathd2<=adate+90

gen xd6=0
replace xd6=1 if deathd2<=adate+30

gen xd7=0
replace xd7=1 if deathd2<=adate+7

gen xd8=0
replace xd8=1 if deathd2<=adate+1

forval x=1/8 {
replace xd`x'=d`x' if year(deathd2)<2009
}


forval x=4/8 {
rename d`x' d`x'_gun
gen d`x'=d`x'_old
}

forval x=1/3 {
rename d`x' d`x'_gun
gen d`x'=d`x'_old
}

replace d1=d1_gun if adate+1095>date("20171231", "YMD")
replace d2=d2_gun if adate+730>date("20171231", "YMD")
replace d3=d3_gun if adate+365>date("20171231", "YMD")

forval x=1/8 {
replace d`x'=. if deathd2<date("20090101", "YMD")
}


*********************************************
* Age at ambulance dispatch
* Birth date is in year-month format
*********************************************

tostring year, gen(yr)
tostring month, gen(mo)
replace mo="0"+mo if month<10
gen str6 mamb=yr+mo

gen adat2=date(mamb, "YM")
gen bdat=date(fmanad, "YM")

gen age=adat2-bdat

replace age=. if age<0

save "C:Data\data2.dta", replace


************************************************
* Start Analysis *******************************
************************************************


*********************************************
* Heterogeneity
* FIGURE III in article
* TABLE B.9 in OL Appendix
*********************************************

cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

foreach x in d1 {

local k=1
est clear

rename p1 var
rename popri inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter popri
rename var p1

rename female var
rename spri inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter spri
rename var female

rename old var
rename alpri inter
qui reghdfe `x' inter private if amb==1 & extra==0, absorb(fe850 year month priout alder) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter alpri
rename var old

rename hvis2 var
rename prihvis inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter prihvis
rename var hvis2

rename multis var
rename primult inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter primult
rename var multis

rename pcvd var
rename pripcvd inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter pripcvd
rename var pcvd

rename presp var
rename pripresp inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter pripresp
rename var presp

rename ppsy var
rename prippsy inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1
rename inter prippsy
rename var ppsy

rename psur var
rename pripsur inter
qui reghdfe `x' inter private var if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum var if e(sample)==1
estadd scalar mean=r(mean)
local k=`k'+1

lab var inter "Private x Variable"
lab var private "Private ambulance"
lab var var "Variable"

* Figure

coefplot (A1 A2 A3 A4 A5 A6 A7 A8 A9), xline(0) keep(inter) aseq swapnames coeflabels(A1 = "Priority 1" A2 = "Male" A3 = "Above age 70" A4="No Hospital visits 2007-2008" A5="Multiple diagnoses  2007-2008" A6 = "CVD  2007-2008" A7 = "Respiratory  2007-2008" A8 ="Psychiatric/drugs  2007-2008" A9="Surgery  2007-2008") graphregion(color(white)) ylabel(,labsize(small)) msymbol(O) xlab(, nogrid)

graph export "input/`x'hetfig.pdf", as(pdf) replace

* Table 

esttab A1 A2 A3 A4 A5 A6 A7 A8 A9 using input/`x'hettab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private inter var) 

esttab A1 A2 A3 A4 A5 A6 A7 A8 A9 using input/`x'hettaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Interaction variable mean" "Observations"))
rename inter pripsur
rename var psur

}


*********************************************
* Grid size fe figure + standard errors
* FIGURE C.1 in OL Appendix
*********************************************
cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

seq fie, f(200) b(50) t(1000)
egen tqq=tag(fie)
replace fie=. if tqq==0

gen beta=.
gen cih=.
gen cil=.
gen enob=.

forval z=200(50)1000 {
qui reghdfe d1 private if amb==1 & extra==0 , absorb(fe`z' year month priout) vce(cluster fe`z')
replace beta=_b[private] if fie==`z'
replace cih=_b[private]+1.96*_se[private] if fie==`z'
replace cil=_b[private]-1.96*_se[private] if fie==`z'
}

lab var fie "Fixed Effects: (Square root) Number of Grids"

tw rarea cih cil fie , astyle(ci) fintensity(30) || line beta fie , lstyle(ci) || scatter beta fie ,  mcolor(black) yline(0) graphregion(color(white)) legend(off) xlabel(200(100)1000) xtitle("Fixed Effects: (Square Root) Number of Grids", margin(0 0 0 3)) msymbol(O) xlab(, nogrid)

graph export ../paper/input/gridfe.pdf, as(pdf) replace

*********************************
* 2 year mortality results
* FIGURE D.3 in OL Appendix
*********************************

cd "C:\Data"

use "C:\data2.dta", clear
set scheme plotplain
seq fie, f(200) b(50) t(1000)
egen tqq=tag(fie)
replace fie=. if tqq==0

gen beta=.
gen cih=.
gen cil=.
gen enob=.

forval z=200(50)1000 {
qui reghdfe d2 private if amb==1 & extra==0 , absorb(fe`z' year month priout) vce(cluster fe`z')
replace beta=_b[private] if fie==`z'
replace cih=_b[private]+1.96*_se[private] if fie==`z'
replace cil=_b[private]-1.96*_se[private] if fie==`z'
}

lab var fie "Fixed Effects: (Square root) Number of Grids"

tw rarea cih cil fie , astyle(ci) fintensity(30) || line beta fie , lstyle(ci) || scatter beta fie ,  mcolor(black) yline(0) graphregion(color(white)) legend(off) xlabel(200(100)1000) xtitle("Fixed Effects: (Square Root) Number of Grids", margin(0 0 0 3)) msymbol(O) xlab(, nogrid)

graph export ../paper/input/gridfe2.pdf, as(pdf) replace

***********************************************
* Randomization inference 
* FIGURE C.2 in OL Appendix
***********************************************
cd "C:\Data"

use "data2.dta", clear

ritest private _b[private], reps(5000) seed(125) strata(fe850) saving("ri.dta", replace) : reghdfe d1 private if amb==1 & extra==0, absorb(fe850 year month priout)

use "ri.dta", clear 

gen beta=abs(_pm_1)
lab var beta "Absolute value of coefficients on Private"
hist beta, width(0.0001) xline(0.0042) graphregion(color(white))
graph export "input\ri_beta.pdf", as(pdf) replace


**************************************
* Post 2011 data only
* TABLE B.5 in OL Appendix
**************************************

cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

foreach x in resp top d1 {

qui reghdfe `x' private if amb==1 & extra==0 & year>2011, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum `x'  if e(sample)==1 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 A3 using input/ltab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 using input/ltaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


***********************************
* Close to Stockholm city sample
* TABLE B.3 in OL Appendix
***********************************

cd "C:\Data"

use "data2.dta", clear

gen dist=0
replace dist=1 if STATION1=="Gustavsberg" | STATION1=="Hallstavik" | STATION1=="Handen" | STATION1=="M�rsta" | STATION1=="Norrt�lje" | STATION1=="Nyn�shamn" | STATION1=="S�dert�lje" | STATION1=="UpplandsV�sby" | STATION1=="Vallentuna" | STATION1=="V�rmd�" | STATION1=="V�ster Haninge" | STATION1=="V�sterhaninge" | STATION1=="�kersberga" | STATION1=="Tyres�"  | STATION1=="Reservambulans"  | STATION1=="Bro"  | STATION1=="Extra"  | STATION1=="Eker�" | STATION1=="Botkyrka"

* Keep city, s�dermalm, h�gersten, farsta, nacka, t�by, jakobsberg, solna, sollentuna, 

local k=1

est clear

foreach x in 8 7 6 5 4 3 2 1  {

qui reghdfe d`x' private if amb==1 & extra==0 & dist==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum d`x'  if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/dmtab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/dmtaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


**********************************************************
* Descriptive statistics
* TABLES B.1 and B.2 in OL Appendix
**********************************************************

cd "C:\Data"

use "data2.dta", clear

* Summary 1 pre-treatment

lab var ambly "Ambulance within 365 d"
lab var alder "Age"
lab var female "Female"
lab var sb "Nursing home resident"
lab var pc_vis "Health center visits"
lab var ph_vis "Hospital visits"
lab var hvis "Any hospital visit"
lab var vcvi "Any health center visit"
lab var dia "Diabetes"
lab var hyp "Hypertension"
lab var heart "Heart failure"
lab var kol "COPD"
lab var angi "Angina pectoris"
lab var opi "Opioid disease"
lab var dep "Depressive disease"
lab var alc "Alcohol disease"
lab var reu "Reumatic"
lab var pain "Cronic pain"
lab var anemi "Anemia"
lab var flim "Cardiac fibrillation"
lab var heart2 "Chronic heart disease"
lab var uvi "Urinary infection"
lab var hyty "Hyperthyroidism"
lab var hylip "Hyperlipidaemia"
lab var stroke "Stroke"
lab var pneu "Pneumonia"
lab var pros "Benign prostate cancer"
lab var back "Back pain"
lab var anx "Anxiety"
lab var pros2 "Malign prostate cancer"
lab var more2 "Multiple units"
lab var pcvd  "Group CVD"
lab var presp "Group Respiratory"
lab var ppsy "Group Psychatric/drugs"
lab var psur "Group Surgery"
lab var multis "Group Diagnoses $>$1"


* Predetermined before transport
local summa1 "ambly alder female sb more2"

* Health 2007-2008
local summa2 "vcvi hvis pc_vis ph_vis dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 pcvd presp ppsy psur multis"

est clear

forval x=1/2 {
qui estpost su `summa`x'' 
est store A`x'
}

esttab A1 using "input/sum1.tex", replace refcat(ambly "\bf{Pre-determined characteristics}", nolabel) cells("mean(fmt(3) label(\bf{Mean})) sd(fmt(3) label(\bf{Std.})) min(fmt(2) label(\bf{Min})) max(fmt(2) label(\bf{Max})) count(fmt(0) label(\bf{Obs}))") label substitute(# \#) fragment nomtitle nonumber noobs 

esttab A2 using "input/sum1.tex", append refcat(vcvi "\addlinespace \bf{Diagnoses}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)


* Summary 2 - outcomes and treatment

cd "C:\Data"

use "data2.dta", clear

est clear

lab var private "Private ambulance"
lab var d1 "Mortality 3 years"
lab var d2 "Mortality 2 years"
lab var d3 "Mortality 1 year"
lab var d4 "Mortality 6 months"
lab var d5 "Mortality 3 months"
lab var d6 "Mortality 1 month"
lab var d7 "Mortality 1 week"
lab var d8 "Mortality 1 day"
lab var resp "Response time (s)"
lab var top "Time to patient (s)"
lab var ptime "Time at patient (s)"
lab var toh "Time to hospital (s)"
lab var fika "Time at hospital (s)"
lab var cond "Assessed severity"
lab var prioin "Priority to hospital"
lab var pqt "Time to ambulance dispatch (s)"
lab var priout "Priority to patient"
lab var pmed "Patient left at home"

lab var ospec "Unspecified/general" 
lab var bcvd "CVD"
lab var bpsy "Alcohol abuse/Psychiatry"
lab var bresp "Respiratory"
lab var bsur "Surgery"
lab var ca "Cardiac arrest"

lab var Aicd2 "Infectious diseases"
lab var Cicd2 "Cancer"
lab var Eicd2 "Endocrinology (Diabetes)"
lab var Ficd2 "Alcohol abuse/Psychiatry"
lab var Gicd2 "Neurology"
lab var Iicd2 "CVD"
lab var Jicd2 "Respiratory diseases"
lab var Kicd2 "Intestinal diseases"
lab var Micd2 "Arthritis/Orthopedics"
lab var Nicd2 "Urology"
lab var Ricd2 "Surgery"
lab var Sicd2 "Minor trauma"
lab var Ticd2 "Poisoning/burn/fracture"
lab var Zicd2 "Without medical condition"


* Treatment
local summa1 "private"

* Controls
local summa2 "priout"

* Time stamps
local summa3 "pqt resp top cond pmed "

local summa4 "ospec bcvd bpsy bresp bsur ca"
* Health outcomes
local summa5 "d1 d2 d3 d4 d5 d6 d7 d8"

* Hospital admissions
local summa6 "Aicd2 Cicd2 Eicd2 Ficd2 Gicd2 Iicd2 Jicd2 Kicd2 Micd2 Nicd2 Ricd2 Sicd2 Ticd2 Zicd2"

est clear

forval x=1/6 {
qui estpost su `summa`x'' 
est store A`x'
}

esttab A1 using "input/sum2.tex", replace refcat(private "\bf{Treatment}", nolabel) cells("mean(fmt(3) label(\bf{Mean})) sd(fmt(3) label(\bf{Std.})) min(fmt(2) label(\bf{Min})) max(fmt(2) label(\bf{Max})) count(fmt(0) label(\bf{Obs}))") label substitute(# \#) fragment nomtitle nonumber noobs 

esttab A2 using "input/sum2.tex", append refcat(priout "\addlinespace \bf{Controls}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)

esttab A3 using "input/sum2.tex", append refcat(pqt "\addlinespace \bf{Outcomes: Ambulance responses}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)

esttab A4 using "input/sum2.tex", append refcat(ospec "\addlinespace \bf{Outcomes: Ambulance Diagnoses}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)

esttab A5 using "input/sum2.tex", append refcat(d1 "\addlinespace \bf{Outcomes: Health}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)

esttab A6 using "input/sum2.tex", append refcat(Aicd2 "\addlinespace \bf{Outcomes: Hospital admission diagnoses}", nolabel) cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") label substitute(# \#) fragment nomtitle nonumber noobs plain collabels(none)


*****************************************************
* Balance of predetermined characteristics
* FIGURE I and FIGURE II in article
*****************************************************

cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

lab var hovs "Hospital visits normalized"
lab var kol "COPD"
lab var angi "Angina pectoris"
lab var opi "Opioid disease"
lab var dep "Depressive disease"
lab var alc "Alcohol disease"
lab var reu "Reumatic"
lab var pain "Cronic pain"
lab var anemi "Anemia"
lab var flim "Cardiac fibrillation"
lab var heart2 "Chronic heart disease"
lab var uvi "Urinary infection"
lab var hyty "Hyperthyroidism"
lab var hylip "Hyperlipidaemia"
lab var stroke "Stroke"
lab var pneu "Pneumonia"
lab var pros "Benign prostate cancer"
lab var back "Back pain"
lab var anx "Anxiety"
lab var pros2 "Malign prostate cancer"
lab var more2 "Multiple units"

local k=1
est clear
foreach x in asd female vcvi hvis hovs sb ambly pvis more2 {
qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
*qui sum `x'
*estadd scalar mean=r(mean)
local k=`k'+1
}

coefplot (A1 A2 A3 A4 A5 A6 A7 A8 A9), xline(0) keep(private) aseq swapnames coeflabels(A1 = " " A2 = "                  " A3 = " " A4 = "                    " A5 = " " A6 =" " A7=" " A8=" " A9=" " A10=" ") graphregion(color(white)) ylabel(,labsize(small)) xscale(range(-0.08(0.02)0.02)) xlabel(-0.08(0.02)0.02) title("Design based") msymbol(O) xlab(, nogrid)

graph save Graph "input\figbal2a.gph", replace


local k=10

foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 pcvd presp ppsy psur multis {
qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
local k=`k'+1
}

coefplot (A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22 A23 A24 A25 A26 A27 A28 A29 A30 A31 A32 A33 A34 A35 A36), xline(0) keep(private) aseq swapnames coeflabels(A9 ="                      " A10=" " A11 = " " A12 = " " A13 = "                      " A14 = " " A15 =" " A16 =" " A17=" " A18=" " A19=" " A20=" " A21=" " A22=" " A23=" " A24=" " A25=" " A26=" " A27=" " A28=" " A29=" " A30="                      " A31=" " A32=" " A33=" " A34=" " A35=" " A36=" ") graphregion(color(white)) ylabel(,labsize(small)) xscale(range(-0.01(0.005)0.015)) xlabel(-0.01(0.005)0.015) title("Design based") msymbol(O) xlab(, nogrid)

graph save Graph "input\figbal2b.gph", replace

********************************************
* Second unconditional: 2 panels
********************************************

* First patient characteristics

local k=1
est clear
foreach x in asd female vcvi hvis hovs sb ambly pvis more2 {
qui reg `x' private if amb==1 & extra==0, cluster(fe850)
est store A`k'
*qui sum `x'
*estadd scalar mean=r(mean)
local k=`k'+1
}

coefplot (A1 A2 A3 A4 A5 A6 A7 A8 A9), xline(0) keep(private) aseq swapnames coeflabels(A1 = `""Age" "(normalized)""' A2 = "Female" A3 = `""Any primary" "care visit""' A4 = `""Any hospital" "visit""' A5 = "Hospital visits" A6 =`""Nursing Home" "Resident""' A7=`""Previously" "ambulance" "transported""' A8=`""Primary care" "visits""' A9=`""Multiple" "units"') graphregion(color(white)) ylabel(,labsize(small))  title("Bivariate (Correlation)") msymbol(O) xlab(, nogrid)

graph save Graph "input\figbal1a.gph", replace

* Then diagnoses for unconditional

local k=10
est clear
foreach x in dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 pcvd presp ppsy psur multis {
qui reg `x' private if amb==1 & extra==0, cluster(fe850)
est store A`k'
*qui sum `x'
*estadd scalar mean=r(mean)
local k=`k'+1
}

coefplot (A10 A11 A12 A13 A14 A15 A16 A17 A18 A19 A20 A21 A22 A23 A24 A25 A26 A27 A28 A29 A30 A31 A32 A33 A34 A35 A36), xline(0) keep(private) aseq swapnames coeflabels(A10 ="Diabetes" A11="Hyper tension" A12 = "Heart failure" A13 = "COPD" A14 = "Angina pectoris" A15 = "Opioid abuse" A16 ="Depression" A17 ="Alcohol abuse" A18="Reumathic" A19="Cronic Pain" A20="Anemia" A21="Cardiac fibrillation" A22="Chronic heart disease" A23="Urinary infection" A24="Hyperthyroidism" A25="Hyperlipidaemia" A26="Stroke" A27="Pneumonia" A28="Benign prostate cancer" A29="Back pain" A30="Anxiety" A31="Malign prostate cancer" A32="CVD" A33="Respiratory" A34="Psychiatric/Drugs" A35="Cancer/pain" A36="Multiple diagnoses") graphregion(color(white)) ylabel(,labsize(small)) title("Bivariate (Correlation)") msymbol(O) xlab(, nogrid)


graph save Graph "input\figbal1b.gph", replace

* Combine figures

graph combine  "input/figbal1a" "input/figbal2a" , ycom cols(2)  graphregion(color(white))
graph export "input\fbal1.pdf", as(pdf) replace
graph export "input\fbal1.png", as(png) replace


graph combine  "input/figbal1b" "input/figbal2b" , ycom cols(2)  graphregion(color(white))
graph export "input\fbal2.pdf", as(pdf) replace
graph export "input\fbal2.png", as(png) replace



****************************************************
* Density distribution in age (important for mortality)
* FIGURE A.2 in OL Appendix
****************************************************
* First, with these narrow FE, there are more cells
* without common support of treatment. These cells
* will end up with age=0 and are more common for
* private that operate rural districts. I condition these
* cells out by taking the mean of private not equal to 1.

cd "C:\Data"

use "data2.dta", clear

egen mean2=mean(private), by(fe850)

drop res
qui areg alder i.year i.month i.priout if amb==1 & extra==0 & mean2!=0 & mean2!=1, absorb(fe850)
predict res if e(sample)==1, res

lab var alder "Patient Age"

twoway kdensity alder if private & res!=., lpattern(dash) || kdensity alder if !private & res!=., legend(label(1 "Private") label(2 "Public")) saving(age1, replace) graphregion(color(white)) title("Unconditional density") xtitle(Age) ytitle("")

twoway kdensity res if private & res!=., lpattern(dash) || kdensity res if !private & res!=., legend(label(1 "Private") label(2 "Public")) saving(age2,  replace) graphregion(color(white)) title("Conditional density") xtitle(Age) ytitle("")

gr combine "age1" "age2", ycom graphregion(color(white))

graph export "input\aged.pdf", as(pdf) replace
graph export "input\aged.png", as(png) replace


*****************************************************
* Summary statistics, as table
* TABLE B.7 in OL Appendix
*****************************************************
cd "C:\Data"

use "data2.dta", clear


lab var pcvd "Group CVD"
lab var presp "Group Respiratory"
lab var ppsy "Group Psychatric"
lab var psur "Group Surgery"
   
qui reg alder  private if amb==1 & extra==0, vce(cluster NyID)
est store A1

qui reghdfe alder  private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A2

lab var private "`:var label alder'"
esttab A1 A2 using input/btab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none)  cells(b(fmt(4)) & se(par fmt(4))) nomtitles nostar parentheses keep(private) 


est clear

foreach x in female vcvi hvis hovs sb ambly pvis more2 dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 pcvd presp ppsy psur multis {
local k=1

qui reg `x' private if amb==1 & extra==0, vce(cluster fe850)
est store A1

qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A2

lab var private "`:var label `x''"
esttab A1 A2 using input/btab.tex, b(4) se(4) append f plain label booktabs noobs collabels(none)  cells(b(fmt(4)) & se(par fmt(4))) nomtitles nostar parentheses keep(private) 

est clear
}

*****************************************************
* Contracted outcomes
* Table II in article
*****************************************************

cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

foreach x in resp top {

qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum `x'  if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 using input/ctab.tex, b(2) se(2) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 using input/ctaba.tex, b(2) se(2) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(0 0) labels("Public outcome mean" "Observations"))


*****************************************************
* Mortality outcomes
* TABLE III in article
*****************************************************
cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

foreach x in 8 7 6 5 4 3 2 1  {

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum d`x'  if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/mtab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/mtaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


******************************************************
* Mortality conditional on transportation status
* TABLE B.8 in OL Appendix
******************************************************


* 1. Dead and left at home


cd "C:\Data"

use "data2.dta", clear


local k=1
est clear

foreach x in 8 7 6 5 4 3 2 1  {

gen dx`x'=d`x'
replace dx`x'=0 if pmed==0
replace dx`x'=. if pmed==.
qui reghdfe dx`x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum dx`x' if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1
}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/m2tab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/m2taba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))

* 2. Dead and tansported

cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

foreach x in 8 7 6 5 4 3 2 1  {

gen dx`x'=d`x'
replace dx`x'=0 if pmed==1
replace dx`x'=. if pmed==.
qui reghdfe dx`x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum dx`x' if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1
}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/m3tab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6 A7 A8 using input/m3taba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


***********************************************************
* Oster selection on unobservables
* Section C.2 in OL Appendix
***********************************************************

cd "C:\Data"

use "data2.dta", clear

* Residualize all variables (does not work well with all the Fixed Effects)

foreach x in d1 private alder female sb ambly vcvi hvis pc_vis ph_vis dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 more2  {

qui reghdfe `x' if amb==1 & extra==0, absorb(fe850 year month priout) res(q_`x')

}

qui reg q_d1 q_private q_ambly q_vcvi q_hvis q_pc_vis q_ph_vis q_alder q_female q_sb q_more2 q_dia q_hyp q_heart q_kol q_angi q_opi q_dep q_alc q_reu q_pain q_anemi q_flim q_heart2 q_uvi q_hyty q_hylip q_stroke q_pneu q_pros q_back q_anx q_pros2 if amb==1 & extra==0, vce(cluster fe850)

* Use 0.247 as RMAX, 1.3*R^2adj (0.19) from controlled regression.

psacalc delta q_private, beta(0) rmax(0.247)

***********************************************************
* All time stamps
* TABLE B.4 in OL Appendix
***********************************************************
cd "C:\Data"

use "data2.dta", clear

local k=1
foreach x in resp top ptime toh fika {
qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum `x' if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1
}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 using input/stampstab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 using input/stampstaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


***************************************************
* Behaviour of ambulance crews
* TABLE IV in article
***************************************************

cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

foreach x in ospec bcvd bpsy bresp bsur ca cond pmed  {

qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum `x'  if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6 A7 A8  using input/behtab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6 A7 A8  using input/behtaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))

****************************************************
* Hospital admissions
* TABLE V in article
****************************************************

cd "C:\Data"

use "data2.dta", clear

local k=1

est clear

foreach x in inl inl_all Iicd2 Ficd2 Jicd2 Kicd2  {

qui reghdfe `x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
qui sum `x'  if amb==1 & extra==0 & private==0
estadd scalar mean=r(mean)
local k=`k'+1

}

lab var private "Private ambulance"
esttab A1 A2 A3 A4 A5 A6  using input/inltab.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 A3 A4 A5 A6  using input/inltaba.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


*******************************************************
* Diagnosis at intake - consequences of actions taken
* FIGURE F.1 in OL Appendix
*******************************************************

cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

local k=1
foreach z in A C E F G I J K M N R S T Z {
qui reghdfe `z'icd2 private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
eststo A`k'
local k=`k'+1
}

coefplot (A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14), xline(0) keep(private) aseq swapnames coeflabels(A1 = "Infectious diseases (A)" A2 = "Cancer (C)" A3 = "Endocrinology (Diabetes) (E)" A4 = "Alcohol abuse/Psychiatry (F)" A5 = "Neurology (G)" A6 ="CVD (I)" A7="Respiratory diseases (J)" A8="GI (K)" A9="Arthritis/Orthopedics (M)" A10="Urology (N)" A11="Unclassified (R)" A12="Minor trauma (S)" A13="Poisoning/burn/fracture (T)" A14="Without medical condition (Z)") graphregion(color(white)) ylabel(,labsize(small))  msymbol(O) xlab(, nogrid)

graph export "input\inlcause.pdf", as(pdf) replace

*******************************************************
* Diagnosis at intake - CVD decomposed
* FIGURE F.2 in OL Appendix
*******************************************************

cd "C:\Data"

use "data2.dta", clear

* Most common intake diagnoses in heart category:
* I214 (9000) AMI
* I219 (5600) AMI
* I209 (4000) angina pectoris
* I489 (6000) fibrillation
* I509 (18000) heart failure
* I634 (4500) stroke
* I639 (10000) stroke 

set scheme plotplain

gen temp2=substr(diag1,1,3)

* AMI
gen iami=0
replace iami=1 if temp2=="I21"

* Angina
gen iap=0
replace iap=1 if temp2=="I20"

* Fibrillation
gen ifib=0
replace ifib=1 if temp2=="I48"

* Heart failure
gen ihf=0
replace ihf=1 if temp2=="I50"

* Stroke
gen istr=0
replace istr=1 if temp2=="I63"

local k=1
est clear
foreach z in iami iap ifib ihf istr {
qui reghdfe `z' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
eststo A`k'
local k=`k'+1
}

coefplot (A1 A2 A3 A4 A5), xline(0) keep(private) aseq swapnames coeflabels(A1 = "AMI" A2 = "Angina Pectoris" A3 = "Fibrillation" A4 = "Heart failure" A5 = "Stroke") graphregion(color(white)) ylabel(,labsize(small))  msymbol(O) xlab(, nogrid)

graph export "input\intake_cvd.pdf", as(pdf) replace


*****************************************************
* Robustness, shows our identification 
* strategy with many twists that could invalidate it
* FIGURE IV in article and FIGURE D.2 in OL Appendix
*****************************************************

cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

local z=850

forval x=1/2 {

local k=1
est clear

qui reg d`x' private if amb==1 & extra==0, cluster(NyID)
est store A`k'
local k=`k'+1


qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' year month) vce(cluster fe`z')
est store A`k'
local k=`k'+1

* Prefered model
qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private  if amb==1 & extra==0 & aba==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1


qui reghdfe d`x' private  if amb==1 & extra==0 & more==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private  if amb==1 & extra==0 & priout<4, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private  if amb==1 & extra==0 & priout<3, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private alder female sb ambly vcvi hvis pc_vis ph_vis dia hyp heart kol angi opi dep alc reu pain anemi flim heart2 uvi hyty hylip stroke pneu pros back anx pros2 more2 if amb==1 & extra==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z'#year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z'#year#month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z'#year#week priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

* Hospital FE
qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' year month priout hosp) vce(cluster fe`z')
est store A`k'
local k=`k'+1

* No transfers from hospitals
qui reghdfe d`x' private  if amb==1 & extra==0 & fhosp==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1

* Patient queue
qui reghdfe d`x' pqt private if amb==1 & extra==0, absorb(fe`z' year month priout) vce(cluster fe`z')
est store A`k'
local k=`k'+1


* Time
qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' priout year#month) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' priout year#week) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' priout year#week#dow) vce(cluster fe`z')
est store A`k'
local k=`k'+1

qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe`z' priout year#week#dow#ctime) vce(cluster fe`z')
est store A`k'
local k=`k'+1


coefplot (A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18), xline(0) keep(private) aseq swapnames coeflabels(A1 = "Correlation (1)" A2 = "Add Grid+Year+Month FE (2)"  A3 = "Add Priority to Patient (3)" A4 = "No Emergency Car Ambulances (4)" A5 = "No Multiple Ambulance Transports (5)" A6 = "No Priority 4 (6)" A7 = "Only Priority 1-2 (7)"  A8 = "Add all controls (8)" A9 = "FE Interacted with Year (9)"  A10 = "FE Interacted with Year and Month (10)" A11 ="FE Interacted with Year and Week (11)"  A12 ="Hospital FE (12)"  A13 = "No hospital transfers (13)" A14="Control for call time (14)"  A15="Year/Month FE (15)" A16="Year/Month/Week FE (16)" A17="Year/Month/Week/Day FE (17)" A18="Year/Month/Week/Day/Hour FE (18)") graphregion(color(white)) msymbol(O) xlab(, nogrid)

graph export "input/reg`z'd`x'.pdf", as(pdf) name("Graph") replace
graph export "input/reg`z'd`x'.png", as(png) name("Graph") replace

}


***************************************************
* Mortality as figure
* FIGURE A.3 in OL Appendix
***************************************************
cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

est clear
local k=1
forval x=1/8 {
qui reghdfe d`x' private if amb==1 & extra==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
local k=`k'+1
}

coefplot (A8 A7 A6 A5 A4 A3 A2 A1), yline(0) keep(private) vertical aseq swapnames coeflabels(A8 = "1 Day" A7 = "1 Week" A6 = "1 Month" A5 = "3 Months" A4 = "6 Months" A3 = "1 Year" A2 = "2 Years" A1 = "3 Years") graphregion(color(white)) ciopts(recast(rcap)) msymbol(O) xlab(, nogrid)
graph export "input/reg_te.pdf", as(pdf) name("Graph") replace
graph export "input/reg_te.png", as(png) name("Graph") replace

*****************************************************
* Mortality effect difference between 
* private and public areas?
* TABLE D.1 in OL Appendix
*****************************************************

cd "C:\Data"

use "data2.dta", clear

local k=1
est clear

qui reghdfe d1 private if amb==1 & extra==0 & apriv==0, absorb(fe850 year month priout) vce(cluster fe850)
est store A1
qui sum d1  if amb==1 & extra==0 & private==0 & apriv==0
estadd scalar mean=r(mean)

qui reghdfe d1 private if amb==1 & extra==0 & apriv==1, absorb(fe850 year month priout) vce(cluster fe850)
est store A2
qui sum d1  if amb==1 & extra==0 & private==0 & apriv==1
estadd scalar mean=r(mean)


lab var private "Private ambulance"
esttab A1 A2 using input/tab_area.tex, b(4) se(4) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep(private) 

esttab A1 A2 using input/tab_areaa.tex, b(3) se(3) replace f plain label booktabs noobs collabels(none) nomtitles parentheses keep("") stats(mean N, fmt(4 0) labels("Public outcome mean" "Observations"))


*************************************************************
* Overlap test. Is the effect different where there is
* balanced common support in private-public?
* FIGURE D.1 in OL Appendix
*************************************************************

cd "C:\Data"

use "data2.dta", clear

set scheme plotplain

egen mean2=mean(private), by(fe850 year)

est clear
local k=1
forval x=0.1(0.1)0.4 {
local z=1-`x'
qui reghdfe d1 private if amb==1 & extra==0 & mean2>`x' & mean2<=`z', absorb(fe850 year month priout) vce(cluster fe850)
est store A`k'
local k=`k'+1
}

coefplot (A1 A2 A3 A4), vertical yline(0) keep(private) aseq swapnames coeflabels(A1 = "0.1-0.9" A2 = "0.2-0.8" A3 = "0.3-0.7" A4 = "0.4-0.6") graphregion(color(white)) legend(off)  mlabposition(1) mlabgap(*2) mlabel(string(@b,"%9.2g")) mlabsize(vsmall) xtitle("Sample Included by Grid Cell Prevalence of Private Ambulances", margin(0 0 0 4)) ytitle(Point Estimate on Private Dummy) msymbol(O) xlab(, nogrid)
graph export "input/overlap.pdf", as(pdf) name("Graph") replace


clear

exit

*****************************************
* END CODE FILE
*****************************************




